Raster data in R

This document is a ‘cheatsheet’ for working with raster datasets in R and introduces the terra package which is the preferred option for handling such data currently. We first load some libraries (including some vector libraries too).

library(terra)
library(sf)
library(tmap)
library(dplyr)
library(tidyr)
library(ggplot2)
tmap_mode("plot")

An example dataset to start

Reading a dataset is pretty simple, especially if it is a simple dataset. I have put a bunch of GeoTIFFs in this zip file, which you should download and unzip to a folder on your computer. Then read the digital elevation model dataset dem.tif into a variable with the rast() function, like this:

dem <- rast("dem.tif")

Confusingly this doesn’t tell you anything about what it just did. So you have to ask:

dem
## class       : SpatRaster 
## dimensions  : 242, 240, 1  (nrow, ncol, nlyr)
## resolution  : 100, 100  (x, y)
## extent      : 1735056, 1759056, 5419610, 5443810  (xmin, xmax, ymin, ymax)
## coord. ref. : NZGD2000 / New Zealand Transverse Mercator 2000 (with axis order normalized for visualization) 
## source      : dem.tif 
## name        :     dem 
## min value   :       0 
## max value   : 518.799

This tells us a few things about the dataset, including the dimension (in numbers of cells), the resolution (in units of the projection), the extent, and projection . It also tells us the minimum and maximum values in the data.

As usual in R we can get a summary of the statistics of the dataset:

summary(dem)
##       dem        
##  Min.   :  0.00  
##  1st Qu.: 82.07  
##  Median :149.17  
##  Mean   :156.96  
##  3rd Qu.:219.22  
##  Max.   :518.80  
##  NA's   :23790

What are all those NA values? Well… best bet is to plot the data to make sense of it

plot(dem)

This shows us right away that a lot of the rectangular area of the raster is taken up by sea, and since it is an elevation model, we might expect the sea areas to be recorded as NA values (which they are).

tmap can consume rasters too!

tmap knows about raster data and can map them, quite happily. The same options such as the number of classes (n), the style of classification (e.g. quantile), and the colour palette to us in the mapping also work for rasters. Here, because you might want to use it, I’ve also shown how to use R’s built-in terrain colour palette.

tm_shape(dem) +
  tm_raster(n = 10, palette = terrain.colors(10)) +
  tm_legend(legend.outside = TRUE)
## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum New Zealand Geodetic Datum 2000 in Proj4
## definition

tmap can also map raster layers in combination with non-raster layers. For example, read in the nz.gpkg file also found in the provided zip file, and add it as an outline around the raster:

welly <- st_read("welly.gpkg")
## Reading layer `welly' from data source 
##   `/home/osullid3/Documents/teaching/Geog315/labs/scratch/rasters/welly.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 79 features and 64 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 1735078 ymin: 5419585 xmax: 1759042 ymax: 5443764
## Projected CRS: NZGD2000 / New Zealand Transverse Mercator 2000 (with axis order normalized for visualization)
tm_shape(dem) + 
  tm_raster() + 
  tm_shape(welly) + 
  tm_borders(col = "blue")
## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum New Zealand Geodetic Datum 2000 in Proj4
## definition

But wait, there’s more than one layer!

There are a bunch of raster data sets in this folder. We can read them into a raster stack by making a directory listing of them and supplying that to the rast function

layer_sources <- dir(pattern = "*.tif")
layers <- rast(layer_sources)

Again, you have to ask to get more information

layers
## class       : SpatRaster 
## dimensions  : 242, 240, 11  (nrow, ncol, nlyr)
## resolution  : 100, 100  (x, y)
## extent      : 1735056, 1759056, 5419610, 5443810  (xmin, xmax, ymin, ymax)
## coord. ref. : NZGD2000 / New Zealand Transverse Mercator 2000 (with axis order normalized for visualization) 
## sources     : age.tif  
##               deficit.tif  
##               dem.tif  
##               ... and 8 more source(s)
## names       :        age,    deficit,        dem,        mas,        mat,      r2pet, ... 
## min values  :     0.0000,     0.0000,     0.0000,   138.0000,   100.0000,    24.0000, ... 
## max values  :    2.00000,  133.02213,  518.79901,  143.00000,  131.00000,   41.00000, ...

If we want to know the names of all the layers, we can use the names() function

names(layers)
##  [1] "age"     "deficit" "dem"     "mas"     "mat"     "r2pet"   "rain"   
##  [8] "slope"   "sseas"   "tseas"   "vpd"

And plot will make several plots of all the layers (this might be a bit slow).

plot(layers)

tmap can also deal with multiple layers, although it will apply the same classification to all of them unless you override this behaviour:

tm_shape(layers) + 
  tm_raster() + #title = names(layers)) # this fix no longer needed
  tm_layout(legend.position = c("RIGHT", "BOTTOM")) +
  tm_facets(free.scales = TRUE) # different scales on each layer

There might be some complaining about fitting the legends in on this plot, or about projections (this is to do with recent changes in how projections are handled by the proj tools see this post but you can ignore it for now).

To examine a single layer in the stack, just use the $ notation:

tm_shape(layers$mat) +
  tm_raster(palette = "-RdBu", n = 9, title = "Mean annual temp (x10)")
## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum New Zealand Geodetic Datum 2000 in Proj4
## definition

Layers in the raster stack

The layers included now in the layers dataset are:

The layers in the stack were assembled from LENZ data developed by Manaaki Whenua Landcare Research, and available here. Their interpretations are roughly as follows (this information from the disdat package).

Layer name Explanation Units
age Soil age 3 classes (0 to 2): <2000, 2000-postglacial (app. 30,000), and pre-glacial
deficit Annual water deficit mm (higher is greater deficit)
dem Elevation metres
mas Mean annual solar radiation MJ/m2/day
mat Mean annual temperature degrees C * 10
r2pet Average monthly ratio of rainfall and potential evapotranspiration (ratio)
rain annual precipitation mm
slope Slope degrees
sseas Solar radiation seasonality dimensionless
tseas Temperature seasonality degrees C
vpd Mean October vapor pressure deficit at 9 AM kPa

Raster calculation

terra will also allow you to do operations between layers, to perform raster calculations. For example, we can easily do this (nonsensical) calculation

layers$mas + layers$mat - layers$dem * layers$deficit
## class       : SpatRaster 
## dimensions  : 242, 240, 1  (nrow, ncol, nlyr)
## resolution  : 100, 100  (x, y)
## extent      : 1735056, 1759056, 5419610, 5443810  (xmin, xmax, ymin, ymax)
## coord. ref. : NZGD2000 / New Zealand Transverse Mercator 2000 (with axis order normalized for visualization) 
## source      : memory 
## name        :       mas 
## min value   : -14228.81 
## max value   :       274

This makes a new layer by mathematically combining the values in the original layers. When combined in this way raster layers must be aligned with one another and in the same projected coordinate system, and the values at the same locations in each layer contribute to the result at the corresponding location in the result.

This means that we can use raster calculations of this type to perform overlay as described in the next section.

Raster overlay

If we recode each layer by whether or not it is greater than some threshold, then we can logically combine them to find regions that satisfy a desired combination of criteria. Say we wanted places below 250m elevation with slope under 5 degrees:

below_250m <- layers$dem < 250
flat <- layers$slope < 5
result <- below_250m & flat

As usual make some maps

m1 <- tm_shape(below_250m) + 
  tm_raster(palette = "Set1", style = "cat", title = "dem<250m") +
  tm_shape(welly) + 
  tm_borders(col = "white", lwd = 0.2)
m2 <- tm_shape(flat) + 
  tm_raster(palette = "Set1", style = "cat", title = "slope<5deg") + 
  tm_shape(welly) + 
  tm_borders(col = "white", lwd = 0.2)
m3 <- tm_shape(result) + 
  tm_raster(palette = "Set1", style = "cat") + 
  tm_shape(welly) + 
  tm_borders(col = "white", lwd = 0.2)
tmap_arrange(m1, m2, m3, nrow = 1)

There is nothing to stop us performing all these operations in a single step, once we are used to how it works:

r <- (layers$dem < 250) & (layers$slope < 5)

This makes it relatively easy to flexibly perform complex suitability evaluations using terra raster calculation operations.

Focal calculations

The simple calculations in the previous sections are known as local because they apply locally across the values at the same locations in each layer. Two other kinds of calculation are possible focal and zonal.

We can apply functions in focal windows within a single layer and perform calculations based on the values within that window centred on each location across a raster. For example, we might want to determine the max value of raster cells in a 3x3 window.

max_elev <- layers$dem %>%
  focal(fun = max)
plot(max_elev)

A 3x3 window is the default. If you want a larger window specify the size with the w parameter (which must be an odd number)

max_elev <- layers$dem %>%
  focal(w = 11, fun = max)
plot(max_elev)

Notice how a larger window produces a more ‘smoothed’ or generalised result.

Again, we can flexibly combine such operations. For example a measure of relief (or roughness) is the local difference between the maximum and minimum values of that surface

relief <- focal(layers$dem, fun = max) - focal(layers$dem, fun = min)
names(relief) <- "relief"
tm_shape(relief) +
  tm_raster(style = "cont", palette = "Reds", alpha = 0.7)
## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum New Zealand Geodetic Datum 2000 in Proj4
## definition

Zonal calculations

You can also perform calculations based on a set of zones. If you want the results in a raster, then both inputs must be rasters. That means we have to do quite a lot of fiddling around with the data. First we make a zones dataset from some polygons, converting first to terra‘s SpatVector format, then using rasterize to select one attribute from the zones as an identifier for zones. The raster we provide to the rasterize function is a ’template’ for the number of cells and their resolution.

welly_zones <- welly %>%
  as("SpatVector") %>%
  rasterize(relief, field = "ID")
tm_shape(welly_zones) + 
  tm_raster(n = 30, palette = "Set1", legend.show = FALSE) + 
  tm_shape(welly) +
  tm_borders()
## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum New Zealand Geodetic Datum 2000 in Proj4
## definition

We can’t easily map the zones, because there are a large number of unique values.

Once we have zones we can perform a zonal calculation, which will calculate the result of applying a function across all the raster cells in each zone. We may need to change the NA values in the input raster layer to zeros (or some other defined value) so that any missing NA values in a zone don’t prevent calculation of a result for the zone they are in.

zonal_relief <- relief %>%
  replace_na(0) %>% 
  zonal(welly_zones, fun = mean, as.raster = TRUE) %>%
  mask(relief) ## 

The result of this operation is a raster where the value in each cell is the mean relief of all the cells in the zone it is in.

tm_shape(zonal_relief) + 
  tm_raster(n = 10, palette = "RdYlBu", title = "Zonal mean relief") + 
  tm_shape(welly) + 
  tm_borders()
## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum New Zealand Geodetic Datum 2000 in Proj4
## definition

In many situations, it is more appropriate to output the results from zonal calculations as a table of values which can be attached to the original vector data.

zonal_relief_table <- relief %>%
  replace_na(0) %>% 
  zonal(welly_zones, fun = mean)
welly_relief <- welly %>%
  merge(zonal_relief_table)
tm_shape(welly_relief) +
  tm_polygons(col = "relief", style = "cont")

Surface analysis

terra provides an array of commonly used surface analysis methods. These are a better option than figuring out how to calculate them yourself using focal functions.

For example, we can calculate slopes from an elevation model using the terrain function

slope <- terrain(dem)
tm_shape(slope) + 
  tm_raster()
## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum New Zealand Geodetic Datum 2000 in Proj4
## definition

The terrain function offers several options, which you can find out more about in the help using ?terrain. A common use is for calculating hillshade, which can be done like the example below

aspect <- dem %>% 
  terrain("aspect", unit = "radians")
slope <- dem %>%
  terrain("slope", unit = "radians")
hillshade <- shade(slope, aspect, angle = 315)
tm_shape(hillshade) + 
  tm_raster(palette = "Greys", style = "cont", alpha = 0.65, legend.show = FALSE)
## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum New Zealand Geodetic Datum 2000 in Proj4
## definition

Combining raster and vector operations

In the examples for zonal calculations, it was necessary to combine raster and vector data and translate between them.

This is a common requirement, and can sometimes be challenging. Here a few important examples are shown. The key to success is realising that terra has its own vector data format SpatVector and terra functions want vector data to be in this format in order to work. To convert sf vector data to SpatVector use as("SpatVector") and to convert from SpatVector to sf use st_as_sf().

Extract values to points

For example, say we have a set of points, which we want values from a raster layer attached to.

First load the points. Here we have some sampled plants. They have to be projected to match the target raster, and we also restrict them only to the study area using st_filter().

plants <- st_read("tradescantia.gpkg") %>%
  st_transform(crs(layers)) %>% 
  st_filter(welly)
## Reading layer `tradescantia' from data source 
##   `/home/osullid3/Documents/teaching/Geog315/labs/scratch/rasters/tradescantia.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 1635 features and 3 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -176.5629 ymin: -46.9931 xmax: 178.4418 ymax: -34.7562
## Geodetic CRS:  WGS 84

Extract the values of a particular raster with this point set, we have to convert it to SpatVector points, then run terra::extract (we have to specify terra:: because there is another function called extract in tidyr).

mean_t <- layers$mat %>%
  terra::extract(as(plants, "SpatVector"))
as_tibble(mean_t)
## # A tibble: 120 × 2
##       ID   mat
##    <dbl> <dbl>
##  1     1  125.
##  2     2  119 
##  3     3  126.
##  4     4  114.
##  5     5  117 
##  6     6  130 
##  7     7  121 
##  8     8  123.
##  9     9  125 
## 10    10  120.
## # … with 110 more rows

The result is just a table. If we want to map this, we have to add it back to the original spatial data:

plants <- plants %>%
  mutate(mean_t = mean_t$mat)
tm_shape(plants) + 
  tm_dots(col = "mean_t")

We can apply terra::extract to a raster stack to get a multi-attribute table all in one go:

layers %>% terra::extract(as(plants, "SpatVector")) 
##      ID age   deficit        dem mas      mat    r2pet     rain     slope sseas
## 1     1   2  62.28255  81.850563 142 125.0223 29.00000 1181.559 24.977749    44
## 2     2   2  26.64497 170.385193 141 119.0000 35.00000 1206.175  7.644967    44
## 3     3   2  48.87741  74.046143 141 125.9484 32.05157 1193.613 23.967861    44
## 4     4   2  14.00000 270.647339 139 113.9835 39.00000 1171.995 17.677706    44
## 5     5   2  25.39118 200.619965 141 117.0000 36.00000 1153.038  6.008887    44
## 6     6   2  98.96744   8.207642 143 130.0000 26.00000 1115.785  1.000000    45
## 7     7   2  49.60237 141.479630 142 121.0000 32.00000 1134.068 14.609330    44
## 8     8   2  41.35058 115.502556 140 122.9810 35.00000 1158.736 15.000000    44
## 9     9   2  45.62231  83.624641 141 125.0000 33.00000 1190.220 13.034999    44
## 10   10   2  26.95165 167.459930 140 119.9517 36.00000 1197.365 15.633889    44
## 11   11   2  44.27446  92.351151 141 124.9688 33.00000 1193.510 24.991325    44
## 12   12   2  45.64238  75.081467 141 125.0000 34.00000 1183.438  9.397834    44
## 13   13   2  35.00000 131.923416 141 122.0000 34.00000 1201.759  2.625707    44
## 14   14   2 109.00000   2.933569 143 130.0000 25.00000 1137.889  5.926316    45
## 15   15   2  68.56382  50.979637 142 127.0000 29.00000 1148.681 12.436181    44
## 16   16   2  66.77875  78.654251 142 125.0000 29.00000 1131.083 15.435949    44
## 17   17   2  14.00681 244.173737 140 115.0068 39.00000 1160.724 17.981869    44
## 18   18   2  30.37409 154.893311 141 120.0000 34.00000 1224.243 14.041884    44
## 19   19   2  52.31923 121.309883 142 122.6391 31.00000 1159.858 23.655590    44
## 20   20   2  13.62304 230.870667 140 115.0265 38.00000 1183.262 13.911610    44
## 21   21   2  42.73308 105.846130 141 123.0000 34.00000 1163.599 12.867723    44
## 22   22   2  31.01176 138.779724 141 121.0000 36.00000 1171.339  5.000000    44
## 23   23   2  53.10150  45.262615 142 127.0508 32.00000 1192.101 21.796999    44
## 24   24   2  32.07914 150.639267 141 121.0000 34.00000 1234.301  5.029380    44
## 25   25   2  37.02601 114.676880 141 123.0000 35.00000 1184.033 11.641457    44
## 26   26   2  41.03801  90.667397 141 124.0000 34.00000 1181.224  9.061888    44
## 27   27   2  42.03058 123.655266 141 122.0306 34.00000 1163.380  7.030577    44
## 28   28   2  44.29372  83.716370 141 125.0000 34.00000 1182.644 15.311479    44
## 29   29   2  77.61778   8.234794 142 129.0000 29.00000 1154.349 13.351653    44
## 30   30   2  42.03058 123.655266 141 122.0306 34.00000 1163.380  7.030577    44
## 31   31   2  65.62665  83.118309 141 123.9973 32.00000 1156.229 15.753992    44
## 32   32   2  79.00000   8.101830 142 130.0000 29.00000 1152.648  1.000000    44
## 33   33   2  53.00000  51.516697 141 127.0000 32.00000 1190.664  9.076740    44
## 34   34   2  49.05710 111.361656 141 123.0000 32.00000 1192.307  8.326986    44
## 35   35   2  70.97032  65.887962 142 126.0000 29.00000 1128.287  7.624279    44
## 36   36   2  73.60202  28.743425 142 128.9852 29.00000 1150.323  8.360874    44
## 37   37   2  34.58276 132.193329 141 121.9512 35.00000 1196.248 14.582765    44
## 38   38   2  27.04408 161.311722 141 119.3645 35.03992 1178.057 27.404453    44
## 39   39 NaN       NaN        NaN NaN      NaN      NaN 1173.061       NaN   NaN
## 40   40   2  60.37127  64.083328 142 126.0000 31.00000 1154.409  8.419346    44
## 41   41   2  53.52137  55.878609 142 126.9808 32.00000 1193.438 20.835752    44
## 42   42   2  83.57523  52.906380 142 127.0000 27.00000 1119.425 12.982502    45
## 43   43   2  74.15147  88.796715 142 124.9903 28.00000 1108.978 14.013801    45
## 44   44   2  21.00751 184.083023 140 118.0000 37.00000 1161.571  5.980017    44
## 45   45   2  14.00000 221.800766 140 116.0000 39.00000 1169.772 10.629755    44
## 46   46   2  41.76005 106.200203 141 123.9693 33.00000 1189.269 20.429573    44
## 47   47   2  57.42181  61.602367 142 127.0000 32.00000 1168.331 10.370629    44
## 48   48   2  52.00000 121.581322 142 122.0000 32.00000 1137.967  6.771970    44
## 49   49   2  83.73991  30.918926 140 128.0462 30.00000 1192.443 18.597193    44
## 50   50   2  28.01720 179.246902 141 118.0000 36.00000 1153.922 10.371836    44
## 51   51   2  26.68111 159.442719 141 120.0000 36.00000 1182.651 14.334340    44
## 52   52   2  47.61535  97.040489 141 124.0000 33.00000 1162.150  8.384654    44
## 53   53   2  29.00000 156.824371 141 120.0000 35.00000 1208.444  4.940795    44
## 54   54   2  33.77411 166.106613 141 119.0000 35.00000 1149.401 19.612944    44
## 55   55   2  48.94762  84.229973 141 124.9676 33.00000 1166.231 10.333567    44
## 56   56   2  45.33219  90.564613 141 124.9517 33.00000 1187.568 18.712646    44
## 57   57   2  36.40050 181.974335 141 118.0074 34.00000 1141.459 14.385720    44
## 58   58   2  45.24791  18.098948 142 128.3530 32.00000 1266.528 25.660004    44
## 59   59   2  60.30459  22.126696 142 128.3633 31.00000 1186.056 16.968849    44
## 60   60   2  34.73252 144.096893 141 121.0000 35.00000 1163.673 16.732523    44
## 61   61   2  35.36098 143.841354 140 120.3669 35.00587 1158.460 17.531569    44
## 62   62   2  12.00000 283.724396 140 112.0177 40.00000 1159.308 18.017662    44
## 63   63   2  28.63475 241.912872 140 115.0027 36.00000 1153.634 13.628904    44
## 64   64   2  38.00000 111.528107 141 123.0000 34.00000 1202.460  3.632200    44
## 65   65   2  19.92366 208.842209 141 117.0000 35.00000 1243.084 10.923658    44
## 66   66   2  16.00000 226.955856 140 115.9853 38.00000 1169.255 21.000000    44
## 67   67   2  38.26402 116.910995 141 122.9640 33.00000 1208.814 18.080193    44
## 68   68   2  53.00000 154.381088 139 120.0000 33.00000 1176.578  3.332412    44
## 69   69   2  47.64957  95.506989 141 124.0000 33.00000 1165.147 11.000000    44
## 70   70   2  34.96188  21.545456 142 128.0000 33.00000 1250.021 20.615139    44
## 71   71   2  17.35966 179.844604 140 118.0000 38.00000 1172.090 13.047443    44
## 72   72   2  30.00000 146.185760 140 120.0000 36.00000 1162.991  2.000000    44
## 73   73   2  34.34344 126.360092 141 122.0000 35.00000 1188.610 16.956320    44
## 74   74   2  38.02332 106.306831 141 123.0233 35.00000 1180.208  6.962867    44
## 75   75   2  54.00000  69.539719 142 126.0000 32.00000 1168.274  7.594716    44
## 76   76   2  26.98126 213.554977 141 117.0000 35.00000 1144.254  7.608913    44
## 77   77   2  82.23924  33.021381 140 127.9518 30.00000 1193.644 16.364378    44
## 78   78   2  37.62527 149.359497 141 120.0293 34.00000 1158.269  6.029291    44
## 79   79   2  53.95367  40.422386 141 127.9537 32.00000 1189.704  8.982679    44
## 80   80   2  66.00000  19.230110 142 129.0000 30.00000 1176.920  1.000000    44
## 81   81   2  17.00000 210.726776 140 117.0000 38.00000 1167.195 10.244825    44
## 82   82   2  10.00000 305.004395 140 112.0000 38.97747 1209.088 23.085899    44
## 83   83   2  71.25963  32.688908 142 128.0432 29.00000 1168.091 18.400412    44
## 84   84   2  16.59008 210.250671 140 117.0000 38.00000 1189.317 15.659211    44
## 85   85   2  18.44054 202.489563 140 117.0441 37.00000 1190.649 11.088174    44
## 86   86   2  48.94530  72.169975 142 126.0000 32.00000 1199.115 13.619845    44
## 87   87   2  36.69690 121.016777 141 121.3685 34.02003 1167.205 20.599169    44
## 88   88   2  63.77258  40.579601 142 128.0000 31.00000 1166.070 13.613713    44
## 89   89   2  35.73909 107.742851 141 123.0000 35.00000 1181.229 15.319310    44
## 90   90   2  29.37322 164.632294 141 120.0000 35.00000 1225.180 12.000000    44
## 91   91   2  44.42871 143.166260 141 121.0090 33.00000 1145.764 10.428706    44
## 92   92   2  66.23499  84.663139 142 125.0174 30.00000 1132.459 10.431952    44
## 93   93   2  20.62150 195.659027 141 118.0000 35.00000 1248.306 11.000000    44
## 94   94   2  16.00000 194.681168 140 117.9908 38.00000 1170.047 10.103127    44
## 95   95   2  42.03058 123.655266 141 122.0306 34.00000 1163.380  7.030577    44
## 96   96   2  40.01698 103.714996 141 124.0000 34.00000 1194.018  8.586416    44
## 97   97   2  65.00981  62.735703 142 126.0000 30.00000 1144.148 11.147452    44
## 98   98   2  49.03423  85.221207 141 125.0000 33.00000 1166.159  8.986865    44
## 99   99   2  79.17804  66.546135 142 126.0000 28.00000 1112.752 10.821962    45
## 100 100   2  27.96902 170.000351 141 119.0000 36.00000 1166.285 12.667546    44
## 101 101   2  21.32281 167.763794 140 119.0000 37.00000 1172.219 15.029562    44
## 102 102   2  54.71941  52.942146 142 126.3812 32.00000 1181.205 21.899353    44
## 103 103   2  33.38015 185.822006 141 118.0000 34.01198 1143.298 13.607867    44
## 104 104   2  28.37443 150.247620 141 120.0000 36.00000 1171.247  4.374433    44
## 105 105   2  30.33723 154.650818 141 120.0000 35.00000 1181.551 17.973600    44
## 106 106   2  50.09616  61.259510 141 126.0481 33.00000 1191.023 18.797804    44
## 107 107   2  39.42589 115.959381 141 123.0000 34.00000 1187.744 11.260800    44
## 108 108   2  21.90289 135.228958 141 121.9029 34.00000 1252.930 15.555657    44
## 109 109   2  14.01309 252.392471 140 115.0000 39.00000 1182.723 24.097075    44
## 110 110   2  44.95111  82.029579 141 125.0000 33.00000 1193.167 12.902219    44
## 111 111   2  14.00000 272.352417 139 114.0000 39.00000 1172.712 15.649142    44
## 112 112   2  86.42541  59.828053 142 127.0000 27.00000 1102.239 12.013014    45
## 113 113   2  37.62217 136.897842 141 121.0000 35.00000 1159.452  5.746593    44
## 114 114   2  54.73362 130.293320 141 122.0014 33.00000 1153.900 27.268824    44
## 115 115   2  41.23255  89.895660 142 124.9531 32.00000 1253.917 23.646996    44
## 116 116   2  79.00000  24.149178 142 129.0000 29.00000 1143.306  2.988217    44
## 117 117   2  38.36826 129.463318 141 122.0000 33.00000 1225.491  9.610627    44
## 118 118   2  34.37909 147.659988 141 121.0000 34.00000 1189.834  8.758190    44
## 119 119   2  78.59055  60.026829 142 126.0154 28.00000 1117.547  9.538243    45
## 120 120   2  47.77022  74.713654 141 125.6732 32.95148 1191.203 23.423868    44
##     tseas      vpd
## 1      55 35.00000
## 2      47 29.00000
## 3      57 31.00000
## 4      47 25.00000
## 5      45 28.00000
## 6      60 36.00000
## 7      48 31.00000
## 8      50 28.00000
## 9      55 31.00000
## 10     47 28.00000
## 11     56 31.00000
## 12     56 30.00000
## 13     51 30.00000
## 14     59 37.00000
## 15     58 34.00000
## 16     55 33.00000
## 17     46 26.00000
## 18     49 30.00000
## 19     51 32.04100
## 20     46 27.00000
## 21     53 30.00000
## 22     48 29.00000
## 23     58 32.00000
## 24     49 30.00000
## 25     52 30.00000
## 26     54 30.00000
## 27     52 30.00000
## 28     56 30.00000
## 29     60 34.00000
## 30     52 30.00000
## 31     52 29.00000
## 32     60 34.00000
## 33     58 32.00000
## 34     53 30.00000
## 35     56 34.00000
## 36     59 34.00000
## 37     51 29.00000
## 38     47 29.00000
## 39    NaN      NaN
## 40     57 32.00000
## 41     58 32.00000
## 42     57 36.00000
## 43     54 33.00000
## 44     45 28.00000
## 45     46 26.00000
## 46     54 31.00000
## 47     57 32.00000
## 48     50 31.00000
## 49     59 29.00000
## 50     45 28.00000
## 51     46 29.00000
## 52     54 31.00000
## 53     48 29.00000
## 54     45 29.00000
## 55     55 31.00000
## 56     55 31.00000
## 57     45 29.00000
## 58     59 34.94745
## 59     60 33.00000
## 60     48 29.00000
## 61     46 28.00000
## 62     46 25.00000
## 63     46 26.00000
## 64     52 29.00000
## 65     46 29.00000
## 66     46 26.00000
## 67     53 31.00000
## 68     49 26.00000
## 69     54 31.00000
## 70     57 34.93906
## 71     45 27.00000
## 72     45 28.00000
## 73     51 29.00000
## 74     53 30.00000
## 75     56 32.00000
## 76     45 28.00000
## 77     59 29.00000
## 78     47 30.00000
## 79     59 32.00000
## 80     60 33.00000
## 81     46 26.00000
## 82     45 25.00000
## 83     59 34.00000
## 84     46 27.00000
## 85     46 27.00000
## 86     57 32.00000
## 87     52 30.00000
## 88     59 33.00000
## 89     52 29.00000
## 90     48 29.00000
## 91     48 30.00000
## 92     55 33.00000
## 93     46 29.00000
## 94     45 27.00000
## 95     52 30.00000
## 96     53 29.00000
## 97     56 32.00000
## 98     55 31.00000
## 99     56 34.00000
## 100    45 29.00000
## 101    45 28.00000
## 102    58 32.00000
## 103    45 29.00000
## 104    46 29.00000
## 105    48 29.00000
## 106    57 31.00000
## 107    53 30.00000
## 108    49 32.00000
## 109    47 26.00000
## 110    56 31.00000
## 111    47 25.00000
## 112    57 34.00000
## 113    47 28.00000
## 114    49 28.00000
## 115    55 33.00000
## 116    60 34.00000
## 117    51 30.00000
## 118    50 30.00000
## 119    56 34.00000
## 120    56 31.00000

This is very useful in statistical model building (which we will get to in a week or two.)

Make a raster layer from vector data

To convert vector data to a raster we convert first to SpatVector then to raster with rasterize

plants %>% as("SpatVector") %>% 
  rasterize(dem) %>%
  tm_shape() + 
  tm_raster(col = "dem", palette = "Greys")
## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum New Zealand Geodetic Datum 2000 in Proj4
## definition